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Simplified non-Navier-Stokes model of turbulent flow 
and its flrst numerical realization in D2 



by Krzysztof Moszynski [^|J 

General remarks on the simplified model 

Let us consider the following integro-differential equation 



J A ^ ^ 



t^ with properly chosen initial and boundary condition (of Dirichlet type). This 

is the equation of the simplified model proposed by Marek Burnat (see [1],[2]) 
with the additional diffusion term —vAxp., where z/ is a (small) positive co- 
efficient. 

1^ Similarly as in papers [1] and [2] , t, x, a and (3 are independent variables, 

^' 

t^ • p{t, X, a) is the a-mass density function, 

g 
L_i • a, /3 are a-velocities, 

^ • the integral term K/^M(-,a, /3) is the so called mixer. 

\Q We assume that 



t e [o,T], r>o, xen, a, ^eA, 

where il and A are rectangles in R^. Here we take 

Q = [0, Li] X [0, L2], A= [D,, Gi] X [A, GJ. 

To get usual Euler quantities it is enough to integrate with respect to the 
variable a: 

• for Euler mass density function 

Q{ti3l) = f^ p{t,x,a)da, 
J A 
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for Euler "impulse density function" (non official terminology!) 



pit,x) = K, ap{t,x,a)da, 
J A 

for Euler impulse of the mass in cu C fi 

imp(t,a;) = / p{t,u)dx, 

for Euler velocity 

p(t,x) J_4^(ip(t,x,a)da 



v{t,x) 



Q{t,x) !ji,p{t-,x,a)da 



We shall define the the function M(-,a, /3) involved in the mixer. Let 

d = \\§\\pit,x,l5) - \\a\\p{t,x,a) 

and 

J ^ for rf > 

^^^^ - j _£;_ for rf < ' 
where || ■ || is the Euclidean norm. Then M is defined by 

r{d)p(t,x,a) if d> 



(2) ^(■'^'^)-lr(rf)p(t,a:,^) if rf < 

It is easy to verify that 

'" M{-,a,l3)dl3da = 0. 
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Equation (1) is equivalent to 

(3) K,[pt + div^(ap) - udiv^Vxp] = i^^ l M{-,a, I3)dl3. 

J A ^ ^ 

hence, integrating both sides of (3) with respect to a we obtain the following 
equation for the Euler quantities: 

(4) git, x) + dwjj){t, x) - udiVa^VxQit, x) = 0, 



Let us integrate both sides of (4) over u G Q. The Gauss Divergence Theorem 
imphes the following relation for the mass m{t,u) = J^^ g{t,x)dx contained 
in u: 

(5) mt{t,uj)+ / np{t,x)dS + u nVxQit,x)dS = 0, 

J dio J dui 

where n is the unit vector external normal to the boundary doo. 

Observe that equation (5) can be read as the the Mass Conservation Law 
for the model considered. 

Any change of the mass in u is possible only as the result of fluxes through 
the boundary doo: 

• of the impulse np(t,x), and/or 

• of the mass nVxQ{t,x). 



First numerical realization of the simplified model in 2D 

Assume that rectangles Q and A are covered by grids Qh and Aah of steps 
hi and ahi i = 1,2 respectively, where 

(6) ^h = {xik,,X2k2), and Aah = {aii^,a2i2) 

where Xj^. = hikj and an. = ahilj, k = {ki,k2), 1= {h, I2), 

h — — h — J ~ ^i 
0<ki< Mi, MRj < L < PRj, i,j = 1, 2. 



I.J ^ Uj 

N 

(7) T, = (t„), tn = nT, n = 0,l,2,---,N. 



The grid for time interval [0, T], T > 0, of the time-step r = ^ is as follows 



For approximation of the function p on the grid (6) (7) let us introduce 
the grid function 

where x ^^^^2 = {hh, h2k2), a i^^i^ = {ahih, 0/12/2)- 

The grid function defined by (8) has to satisfy the following finite difference 

equation 

(9) dul+' + ai<+\,fe,/ + &i<^+i,fe,^ + «2<+fcVi,^ + ^2<+fe+i,z = 

= '^1%; ~ (^l'^ki-l,k2,l_ ^ ^^"^kx+lMl ~ '^2Mfei,fc2-l,i ~ ^2Mfc^^fc2+l,/ + 

"""(Fdi^+i) + F(m")) + din + dir2 
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with 



ahili ufii ahili z//ii 
ai = — Ai , Oi = — Ai 1 , 

4 2 ' 4 2 ' 

0/2,2/2 '^Ai2 , . 0/^2/2 Vp2 

ao = —Ao — : , 02 = — A2 — : , 

2 1 ^ 2 4 2 

(i = 1 + Z/(/ii + /i2), (ii = 1 - v{llx + /i2) 

\ _ ^ \ _ ^ _ ^1 _ ^2 

Ai — — , A2 — 7—, ^1 — 7—, /i2 — TT- 

/il 112 III 112 

The argument m" of the function F is the long (block) vector 

where q = {MRi + Pi?i) (Mi?2 + -P-R2) — 1 (g + 1 is the number of all pairs of 

indexes / = (/i, 12))- Function F(-) is the result of the trapezoidal quadrature 

in D2 of M(-, ■, ■) over A with respect to the variable /3 (see definition of the 

mixer). Note that F is a term non linear with respect to u. 

Dirichlet Boundary conditions are introduced by means of variables diri and 

dir2: 

-ai{DIRL{n) + DIRL{n + l)) ii h = 

din = { -bi{DIRR{n) + DIRR{n + 1)) if h = Mi - 1 , 

else 



-a2{DIRB{n) + DIRB{n + 1)) if fcs 
dir2 = <! -h2\DIRT{n) + DIRT{n + 1)) if k2 

else 
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Here DIRL{-), DIRR{-), DIRB{-), DIRT{-) are Dirichlet conditions at 
the left side, right side at the bottom and on the top side of the rectangle VL 
respectively. Time level is given as argument of DIR ■ (■). 

Equation (9) is simply finite difference approximation of the equation (1) 
on the grid defined in (6) (7). Finite difference approximation is obtained as 
follows: 

• first partial derivatives with respect to variables Xi and X2'- arithmetic 
mean of central finite differences at time levels n and n + 1 

• similarly, second partial derivatives with respect to Xi and X2'- arith- 
metic mean of second (forward-backward) finite differences at time lev- 
els n and n + 1 

• at both time levels n and n + 1 corresponding Dirichlet boundary con- 
ditions are taken into account 

• the mixer term is approximated using the trapezoidal 2D quadrature 
with respect to the variable /3 at time levels n and n + 1; this gives 
corresponding nonlinear terms F. Finally, the arithmetic mean of levels 
n and n -|- 1 is taken. 

Let A and B be the matrices of dimension 

(g + l)MiM2X (g + l)MiM2 

corresponding to the left and right hand side of the linear part of equation 
(9), respectively. We can now write down a compact form of equations (9): 



(10) 



A u 



n+l 



BM"+y[F(M"+l) 



F(m")] + DIR 



where DIR is the sum of all terms introduced by Dirichlet boundary condi- 
tions. Matrices A and B have a block-diagonal structure. 



A 



B 



A 



B 



A 



A 



B 



B 



B. 



each diagonal block corresponds to part of system (9) that depends on a 
fixed / only; hence the block dimension of the matrices A and B is equal to 
q + 1 X q + 1. Matrices A and B of block dimension M2 x M2 have also a 
block structure: 



A 



D B 
A D B 
A D B 



B 



Di -B 

-A Di -B 

-A Di -B 



A D. 
where A, 5, D, Di are of dimension Mi x Mi and 



D 



d hi 
fli d hi 

Oi d hi 



A 



a2 



02 



ai d 



a2 



, Di 



di —hi 
— Oi di —hi 
— fli di 



a2 



B 



-B Di 



-hi 



—ai di 



Solving system (10) 

As is proved in [3], the mixer function M is Lipschitz continuous with respect 
to the variable p. Hence it follows that function F is Lipschitz continuous as 
well. This fact is very important for method of solving the nonlinear system 
(10). In order to solve approximately the system (10), one can apply the 
following simple iteration: 



;ii) 



A u^P+^'^ = B m" 



TK 



[F(u 



(p)) 



F(u")] + DIR, 



where p is the iteration index. Iteration converges for r small enough. At 



any time step n we can put u 



(0) 



u 



and if IIm^p+i) 



u 



ip)\ 



is not too large, 



then we accept u"'~^^ !=s u^P'^^\ At each iteration step it is necessary to solve 
linear system with the (non-symmetric) matrix A: 



(12) 



Ax 



where f is the sum of all known terms in (11), i.e. terms independent of p+1. 
To this end the Richardson iteration was applied: 

(13) Xk+i = Xk + STk, rk = f-Axk 

with optimally chosen relaxation coefficient s. Let us estimate the optimal re- 
laxation coefficient s corresponding to maximum norm || [xi, 0:2, ■ ■ ■ , a;Ar]"^||oo = 
maxj=i 2,...,Ar|a;j|. If x is the solution of linear system (12), then for errors Ck 
and Cfc+i we have 

Cfc+i = (/ - sA)ek. 

Hence 

||6fe-|-l||oO _i II-' "^iilloO ||Cfc||oO. 

Assume 

Oo,o c^o,! 0,0,2 • ai,N 
A= 

_O'N,0 ^Af,! On,2 ■ 0-N,N 

and ai^i > 0. Then it is easy to verify that 

(14) 11/ — sAlloo = max (11 — saj jl + y^ loj J). 

If I^jVi \oi,j\ < Oj.jfor all < i < A^, then J2jyti \oi,j\ + |1 — sai^i\ takes minimal 
value at Sopt = ^- I'^ ^^^ case discussed Cj^j = d > and 

Yl I'^mI ^ 1*^11 + 1*^21 + \bi\ + I&2I, 

hence 

|«i| + \a2\ + \bi\ + I&2 

~d 



/ - SoptAlloo < - — — , ^ = NOR. 



Looking at the coefficients of equation (9), we infer that the sufficient con- 
dition for convergence of the Richardson iteration, NOR < 1 is satisfied if 
proportions of steps in grids (6) (7) are properly chosen, i.e. if Aj = ^ and 
/ij = 1^ for i = 1, 2 are small enough. 

Computational experiment "Collision" 

This experiment was run on the cluster halo2 of the Interdisciplinary Center 
for Mathematical and Computational Modeling of the University of Warsaw. 



Experiment have to be considered as fully "virtual", because coefficients of 
the model have been taken more or less arbitrarily. Till now we had no 
possibility to confront our computational experiments with reality. In such 
situation, the results obtained have only a qualitative character. 
Description of the experiment. 

Four streams of mass enter into the "empty" rectangle Q through its four 
sides. The entering streams are modeled by the Dirichlet boundary conditions 
for the function p, defined on each of four sides of Q by positive functions 
with graphs of triangular shapes and of height increasing in time. Whole 
experiment contains 1000 time steps. After 120 time steps the four 
streams meet in the center of Q and the ffist period of stagnation begins. 
But the mass is always entering into Q and some parts of the streams start 
to go back, rubbing against parts of streams moving in opposite directions. 
After 520 time steps the ffist period of stagnation ends, and ffist eddies 
appear. This period of turbulence ends after 610 time steps, and in this 
moment the second period of stagnation begin, which ends after 870 time 
steps. At this point the second period of turbulence starts with new eddies 
appearing. 

Figures below give certain more interesting stages of this process. The field 
of unit Euler velocity vectors on the central part of Q can be seen. In the 
preprint version of this paper (see < www .mimuw .edu.pl / preprints >) one 
can find more illustrations concerning the experiment " Collision" . 
About the program 

The program that realizes the algorithm described above, was build and 
can be run on the cluster halo2 of the Interdisciplinary Center for Mathe- 
matical and Computational Modeling of the University of Warsaw. Admis- 
sible dimension of the problem depends on the number of processors used. 
Coefficients, initial and boundary conditions, as well as the the number of 
processors can be chosen by the user. 
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